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Abstract 

Background: Boolean network (BN) is a mathematical model for genetic network and control of genetic networks 
has become an important issue owing to their potential application in the field of drug discovery and treatment of 
intractable diseases. Early researches have focused primarily on the analysis of attractor control for a randomly 
generated BN. However, one may also consider how anti-cancer drugs act in both normal and cancer cells. Thus, 
the development of controls for multiple BNs is an important and interesting challenge. 

Results: In this article, we formulate three novel problems about attractor control for two BNs (i.e., normal cell and 
cancer cell). The first is about finding a control that can significantly damage cancer cells but has a limited damage 
to normal cells. The second is about finding a control for normal cells with a guaranteed damaging effect on 
cancer cells. Finally, we formulate a definition for finding a control for cancer cells with limited damaging effect on 
normal cells. We propose integer programming-based methods for solving these problems in a unified manner, 
and we conduct computational experiments to illustrate the efficiency and the effectiveness of our method for our 
multiple-BN control problems. 

Conclusions: We present three novel control problems for multiple BNs that are realistic control models for gene 
regulation networks and adopt an integer programming approach to address these problems. Experimental results 
indicate that our proposed method is useful and effective for moderate size BNs. 




Background 

In the post-genome era, we observe rapid development in 
systems biology, a field focusing on interactions among 
the components of biological systems. Gene regulatory net- 
works (genetic networks, in short) have been proposed to 
better understand the interaction of various kinds of 
genes, proteins and molecules. Many formalisms have 
been developed as models of genetic regulation processes, 
in particular, Boolean network (BN) has received substan- 
tial attention owing to its capacity to capture the switching 
behavior of genetic processes. Furthermore, its dynamic 
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process is rich and complex and can provide insight into 
the global behavior of large genetic regulatory networks. 
A BN is a simple mathematical network: each gene (node) 
takes either 1 (active) or 0 (inactive), and the state of a 
gene is regulated by several genes called its input genes via 
its Boolean functions. It is to be noted that the states of 
genes can be updated synchronously and asynchronously. 
Thus synchronous BNs and asynchronous BNs have been 
proposed to model the two different behaviors. In this 
paper, synchronous BNs are considered since existing 
control studies of BNs are based on synchronous BNs, and 
detection of singleton attractor is equivalent for both mod- 
els. Though synchronous BNs may be too simple com- 
pared with asynchronous BNs, they are effective to model 
and analyze some properties of real gene networks. 
Indeed, they have been used to analyze D. melanogaster 
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embryo development, and the robustness of the genetic 
networks of S. cerevisiae, E. coli, B. subtilis, D. melanoga- 
ster and A. thaliana [1]. 

Many studies have been done on the analysis of 
attractors for randomly given BNs [2,3]- The main rea- 
son is that different attractors can be regarded as differ- 
ent cell types [2]. In [4-7], several approaches have been 
developed to efficiently identify and/or enumerate 
attractors in BNs. It was reported that detection of a 
singleton attractor (i.e., a fixed state) is an NP-hard pro- 
blem [8]. Devloo et al. studied a method by transforma- 
tion to a constraint satisfaction problem [4], while Irons 
proposed another method based on the use of small 
subnetworks [6]. Garg et al. [5] proposed a method 
based on Binary Decision Diagrams (BDDs). However, 
the average-case or worst-case complexity was not ana- 
lyzed theoretically in these studies. Zhang et al. [7] 
developed algorithms to enumerate singleton and small 
attractors, and also studied the time complexities of 
average cases for these algorithms. Furthermore, Akutsu 
et al. [9] developed algorithms with guaranteed worst 
case time complexity for the singleton attractor detec- 
tion for BNs with limited Boolean functions. In addition, 
Datta et al. [10,11] proposed methods for finding a con- 
trol sequence for Probabilistic BNs (PBNs), a probabilis- 
tic extension of BNs. They defined the control problem 
with minimization of the sum of the total control cost 
and the cost of the final state. The cost of control is the 
cost to apply control inputs in some specified states, 
with the higher terminal costs corresponding to undesir- 
able states. Their method is also applicable to finding 
actions of control for BNs, since BNs are special versions 
of PBNs. But, it is necessary for all these approaches to 
deal with 2" x 2" matrices, which makes it difficult to 
apply them to large-scale BNs. Consequently, Chen et al. 
[13] and Kobayashi and Hiraishi [14] proposed integer 
programming based approaches to variants of the PBN 
control problem. Akutsu et al. [15] also proposed an inte- 
ger programming-based approach for constructing the 
attractor control problem for moderate size BNs. The 
integer linear programming-based (ILP) approach is effec- 
tive to determine the optimal solutions subject to a series 
of constraints. Though [15] have shown that attractor 
control is Ej-hard, tne ILP-based method can be applied 
to medium-size BNs. All of these works focused their 
approaches on a single randomly generated BN (i.e., one 
cell type), leaving the analysis of multiple BNs (various 
cell types) unaddressed. However, in real situations, there 
are different kinds of cell types, so it is more realistic to 
perform attractor analysis for multiple BNs. Thus, the 
development of attractor control problems for multiple 
BNs is an important and interesting challenge. Here, we 
propose an integer linear programming approach for 



constructing novel attractor controls for multiple BNs, 
and we provide three novel formulations of the attractor 
control problem to model various realistic cases. We 
consider simultaneous attractor control for multiple BNs 
and, specifically, focus on analyzing attractor control for 
two BNs that each corresponds to normal cells or cancer 
cells. Our objective, motivated by the fact that anti- 
cancer drugs act in both normal and cancer cells, is to 
find the same control (i.e., letting some genes be always 
active (1) and some genes always inactive (0)) for both 
networks. We aim at finding a control that can damage 
cancer cells significantly, while causing only limited 
damage to normal cells. 

As another variant of attractor control, we find single- 
ton attractors for normal cells with a guaranteed dama- 
ging level for cancer cells. In other words, we try to 
investigate whether there exists a control set that can 
damage cancer cells, while at the same time looking for 
a singleton attractor for normal cells under this control. 

It is also of interest to consider finding a singleton 
attractor for cancer cells with limited damage to normal 
cells, to determine if there exists a control set ensuring 
no damage to normal cells, and under such a control 
set, to search for a singleton attractor for cancer cells. 

To utilize ILP, all instances of the original problem are 
converted into ILP, to be able to apply an existing solver. 
These problems are transformed into ILP in a similar 
and systematic way to be shown later. The experimental 
results, using artificially generated BNs and realistic BNs, 
have demonstrated both the efficiency and the effective- 
ness of our proposed method. 

Problem formulations 

In this section, we first give a brief introduction to BNs. 
Since the attractor control problem is based on the 
attractor detection problem (ATTRACTOR DETEC- 
TION), we begin with a brief introduction to that pro- 
blem. We then define three variants of the attractor 
control problem: simultaneous attractor control for two 
BNs (cancer cell vs normal cell), attractor control for 
normal cells under the assumption of significant dama- 
ging cancer cells, attractor control for cancer cells under 
the assumption that normal cells are not damaged. 

Boolean networks 

A Boolean Network (BN) G(V, F ) is a set of vertices 
V = {v 1( v 2 , v„] and a list of Boolean functions 
F = {/i,/ 2 , ...,/„} where/;- : {0, 1}" -> {0, 1}. Define v,(t) 
to be the state (0 or 1) of vertex v, at time t. The rules 
for the regulatory interactions among the genes are then 
expressed by 

Vi{t+ 1) =/i(v(t)), i = 1,2, ••• , n 
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where 

v(0 = [vi(t), ■■■ , v n {i)\ 

is called the Gene Activity Profile (GAP), x A y, x V y, 
x ® y, x is used to represent "AND" of x and j, "OR" of 
x and j, exclusive OR of x and y, and "NOT" of x, 
respectively. We denote IN (v,) as the set of relevant 
input nodes v n , v ik to v,-. The number of relevant 
nodes incoming to v t is called as indegree of v,. K is 
used to represent the maximum indegree of a BN. 

The following is an example of a BN. 

vi(t+ 1) = v 3 (t), 
y 2 (t+ 1) =M0, 
y 3 (t + 1) = fi(f) a t> 2 (f)- 

Each gene v, is updated by a regulatory function/. The 
corresponding truth table of a BN is given in Table 1. 
The dynamics of a BN and state transition diagram are 
shown in Figure 1. For example, the second row of the 
table means that if the state of BN is [0, 0, 1] at time t, 
then the state at time t + 1 will be [1, 1, 0]. Similarly, 
the arc from 001 to 110 signifies that if the state of BN is 
[0, 0, 1] at time t, then the state will be [1, 1, 0] at time 
t + 1. It is easily seen that a BN with n nodes has in total 
2" possible states. Thus, the state transition table and the 
state transition diagram have 2" rows and 2" vertices 
respectively. 

Detection of an attractor 

In a BN, the GAP v(t + 1) at time t + 1 is determined 
by the GAP v(t) at time t. When an initial GAP v(0) is 
given, a BN will finally converge to a set of global states 
(i.e., a directed cycle in the diagram of state transition). 
We call this set an attractor. If an attractor has only 
one state (i.e, v = f(v)), we call it singleton attractor, 
which is corresponding to a fixed point. Otherwise, if it 
consists of p distinct states, i.e., 

v(p)=f(f(.-f(v(0))-)=v(0), 



Table 1 The truth table 



State 


Vi(t) 


v 2 (t) 


v 3 (t) 


vAt + D 


v 2 [t + 1) 


v 3 (t + 1) 


1 


0 


0 


0 


0 


1 


0 


2 


0 


0 


1 


1 


1 


0 


3 


0 


1 


0 


0 


1 


0 


4 


0 


1 


1 


1 


1 


0 


5 


1 


0 


0 


0 


0 


l 


6 


1 


0 


1 


1 


0 


1 


7 


1 


1 


0 


0 


0 


0 


8 


1 


1 


1 


1 


0 


0 



it is called a cyclic attractor of period p. We can see 
from Figure 1 that 010 and 101 are two singleton 
attractors. 

Here, the problem of attractor detection is defined to 
find an attractor for a given BN. It should be noted that 
finding attractors with large periods is not easy. Thus 
we study attractor detection in which upper bound of 
the period is limited by some given threshold pmax. We 
define this problem as follows. Here we consider the 
case of p max = 1 (i.e., the singleton attractor detection 
problem). 

Definition 1: [ATTRACTOR DETECTION] 

Instance: a BN and p max , the maximum period. 

Problem: Output an attractor whose period is at most 
Pmax- If suc h an attractor does not exist, "NULL" should 
be the output. 

For a BN of Table 1, ATTRACTOR DETECTION 
should output either 010 or 101. Since we only consider 
the case of p max = 1 hereafter, we also use v, to denote 
the (steady) state of v h 

Simultaneous attractor control for multiple Boolean 
networks 

A control problem for a single BN was proposed by 
Akutsu et al. [15]. We consider whether there exists a con- 
trol for multiple BNs. We consider two BNs, Ni(V, F{) and 
N 2 (V, F 2 ), representing a cancer cell and a normal cell, 
respectively. Here, V is a set of genes (nodes) while F 1 and 
F 2 are sets of regulation functions (i.e., sets of Boolean 
functions along with input genes) for A/i and N 2 , respec- 
tively. We try to find the same control (i.e., some genes 
always active (1) and some genes always inactive (0)) for 
both networks, considering that anti-cancer drugs act in 
both normal and cancer cells. Therefore, we set our end- 
point as finding a control that causes limited damage to 
normal cells but significant damage to cancer cells. Sup- 
pose control nodes are chosen as Vh • ••• • Vi m and these 
nodes are assigned by Boolean values of b^, ■■■ , b\ m . 
Then, we call v as a singleton attractor if it satisfies the fol- 
lowing conditions (denoted by CONDI): 

(i) v t = bi, if i e {z'i, i m }, 

(ii) V; =fi(y); otherwise. 

Furthermore, in order to evaluate how appropriate 
each attractor state is, we need a score function h(v) for 
cancer cells and g(v) for normal cells from {0, 1}" to the 
set of real numbers. For simplicity, we assume these 
functions are given in linear form as follows: 

h(v) = ^Pa; ■ (1 — u>i) ■ f, and g(v) = ^/i, ■ (1 — wi) ■ v, 

i i 

where w t = 1 (resp., w, = 0) holds if v, is selected 
(resp., not selected) as a control node. This means that 
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we do not take the scores of the selected control nodes 
into account for the score functions (h(v) and g(v)). 
Here, a t and /3, are real constants. For example, if a, = 1 
for all i, then v(t) = (1, 1, 1) is the state with the 
highest score (i.e., the desired state for the cancer cells). 
Similarly, if & = -1 for all i, then v(f) = (0, 0, 0) is the 
state with the highest score (i.e., the desired state for the 
normal cells). We define 

G(v) = h(v)+g(v) 

as the final score function for these two networks. 
Since the singleton attractors may not be uniquely 
determined, considering the worst case, the minimum 
score of singleton attractors is necessary to be maxi- 
mized. Because it is difficult to give a direct formulation 
of this problem, the problem of simultaneous control of 
attractors is defined with 0, which is a threshold of the 
minimum score, as follows. 

Definition 2 : [SIMULTANEOUS ATTRACTOR 
CONTROL] (SAC) 

Instance: 2 BNs, the number of control nodes m, a 
score function G, and a threshold 8. 

Problem: Find m control nodes and 0-1 assignment 
for them where the minimum score of singleton attrac- 
tors is greater than 9. If no such nodes exist, then 
"NULL" is the output. 

We give an example about attractor control for a 
given BN, say, the BN described in Table 1. We assume 
ai = 1, a 2 = 2, a 3 = 3 and m - 1. If v x is a control node 
and 1 is assigned to this node, there exists one singleton 
attractor 101 with score 3. If 0 is assigned to v 2 , two sin- 
gleton attractors 000 and 101 exist. Note that the scores 
of 000 and 101 are 0 and 4, respectively, so the mini- 
mum score is 0. Then, if 0 is assigned to v 3 , 010 is a sin- 
gleton attractor with score 2. Though checking the 



Table 2 Results on Simultaneous Attractor Control 



n/m 


100/10 


120/12 


140/14 


160/16 


time(sec) 


3.41 


7.03 


9.35 


8.36 


#pos/#rep 


5/3.6 


3/4.67 


3/1 


4/2 


#neg/#rep 


3/1 


3/1 


2/1 


3/1 



other three cases is necessary, we can conclude that the 
first case (assigning 1 to v x ) gives the solution. 

Attractor control for normal cells under the assumption 
of significant damage to cancer cells 

We consider another attractor control variant for multi- 
ple networks. We try to investigate whether there exists 
a control that can guarantee significant damage to can- 
cer cells (AAJ and a singleton attractor for normal cells 
(A/2). To ensure that the control can cause significant 
damage to the cancer cells, we introduce a threshold fx 
for A/i to be given later. The score function G(v) 
becomes g(v). It is possible that there exist multiple sin- 
gleton attractors for N 2 , so we maximize the minimum 
score of the singleton attractors and give the threshold 
6 1 for the minimum score. The problem is formulated 
as follows: 

Definition 3: [Attractor Control under Damaging 
Cancer cells substantially] (ACDC) 

Instance: 2 BNs, a score function G, the number of 
control node m, and thresholds 9 1 and £ x , 

Problem: Find 0-1 assignment to m control nodes 
where the minimum score of singleton attractors of N 2 
is no less than 9 lt and the score of the singleton attrac- 
tor for Ni is greater than respectively. If such nodes 
do not exist, then "NULL" is the output. 

Attractor control for cancer cells under the assumption of 
limited damage to normal cells 

We also investigate whether there exists a control that 
ensures limited damage to the normal cells (N 2 ), and 
whether there still exists a singleton attractor for cancer 
cells (N x ). We give a threshold f 2 for the normal cells 
that guarantees keeping the normal cells undamaged and 
convert the score function into G(v)=/z(v). In order to 
obtain the unique singleton attractor for N lt we consider 
maximizing the minimum score function integrating the 
threshold 0 2 for the minimum score of N x . 

Definition 4: [Attractor Control under Keeping Nor- 
mal cells undamaged] (ACKN) 

Instance: 2 BNs, a score function G, the number of 
control node m, and thresholds 0 2 and f 2 , 
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Problem: Find m nodes and a 0-1 assignment of these 
control nodes for which the minimum score of singleton 
attractors for Ni is no less than 6 2 and the score of sin- 
gleton attractors for N 2 is greater than f 2 , respectively. If 
such nodes do not exist, then "NULL" is the output. 

Methods 

Integer programming, in particular, integer linear pro- 
gramming (ILP) is to maximize (or minimize) a linear 
objective function with linear constraints (i.e., linear 
inequalities and linear equations) with all the variables 
taking integer values. From here on, either 0 or 1 is 
assigned to each variable, representing the Boolean 
values. Furthermore, we introduce x t and si to represent 
a 0-1 variable that corresponds to v, for N 1 and N 2 , 
respectively. 

ILP for ATTRACTOR DETECTION 

Here we review how ATTRACTOR DETECTION is for- 
malized as ILP in [15]. Define S b (x) by 



x if b = 1 
x otherwise. 



Then if Boolean function has k inputs, we represent 
it by 

/;(*<!< •" ' **) = v fi(K- ■■■ ' K) A M*fi) A ■ ■ ■ A M***)- 

Since Boolean constraints cannot be used in ILP 
directly, we define Ty(x) as 



x if b = 1 
1 — x otherwise. 



We introduce x i\ -\ in order to represent whether each 
0-1 assignment for {b^, ■■■ , fr;J with fiib^, ■■■ , bi b ) = 1 
is satisfied. If /;(fv • • • , h k ) = 1, we give constraints 

x ihi ... K > (j2 j=1 ,..., h r ^ x 0) ~ ( fe - !)' 

Xi, bil ... K < £E J=1 ,.., fc V X ^' 

in which the former constraint forces x ib il -\ to be 1 if 
Sbiix^) A • • • A Sh k (Xi k ) is satisfied, and the latter forces 
x i\ -\ to be 0 if it is not satisfied. If fi[bi u ■■■ , b, ij = 0, 
a constraint x i,b h -\ = 0 is simply added. These con- 
straints guarantee that 

x i,b h ~\ = 1 

if and only if 

fiiPin ■■■ ' K) A K {x h ) A • • • A S h {x ik ) = 1. 



Finally, for every x it we add the following constraints 

Xi < x i,\-\ 

^■■■^€{0,1}* 

and 

The above constraints are included so as to ensure 
that Xi =f{xi 1 , ■ ■ ■ , X;J holds for each Thus any fea- 
sible solution corresponds to a singleton attractor. A 
simple example is given to illustrate this method. 
Assume x 1 is determined by X\ =/i(x2, X3) = x 2 vx 3 . 
Then fi(x 2 , x 3 ) can be represented as 

/i(x 2 , x 3 ) = (fi(0, 0) Ax 2 Ax 3 ) v (fi(0, 1) Ax 2 Ax 3 ) 

(fl(l,0) AX 2 AX3) V (/i(l, 1) AX 2 AI3) 
= (x 2 A X3) V (x 2 A X3) V (x 2 A X3). 

Thus, this Boolean function can be converted into the 
following inequalities: 

(1 — X 2 ) + (1 — X3) — 1 = — x 2 — x 3 + 1, 
j(l — X 2 + 1 — X3), 

0, 

X 2 + (l — X3) — 1=X 2 — X3, 
j(x 2 + 1 — X3), 
X 2 + X3 — 1, 
\{X 2 +X3), 

^i,oo + *i,oi +^1,10 +Xi,n, 

\{X\fiQ + X\fi\ + Xi,io + *i,u)- 

It is easily seen from this example that the transfor- 
mation is correct. Integrating all the constraints, we can 
formulate the integer programming for singleton attrac- 
tor detection as below. 

maximize x 2 

subject to 

for all i e [1 • ■ n] and 6,, - - • 6,, e |0, l) k , such that /,(!>,„ • • ■ , b it ) = 1, 
^,b, v . = 0, for all i e [1 - ■ • n] and t»; t e 10, l) fe , such that /,(&(,, ■ ■ ■ , &ij = 0, 



*1,00 


> 


*1,00 


< 


Xl,01 




*1,10 


> 


Xl,W 


< 


Xi,n 


> 


Xi,n 


< 


Xi 


< 


Xi 


> 



b.^.-bi €{0,1}' 



^i,b, v ..bi , for all i e [1 • • - n], 



xi e (0,1], forallis [1-n] 



e (0, 11, for all i e [1- n] and b h ■■■li i e (0, 



We note that ATTRACTOR DETECTION is a decision 
problem, not an optimization problem, so we do not need 
an objective function, but we will define an objective func- 
tion in order to apply ILP. Here, we give the objective 
function simply as "Maximize x 2 ". We can extend the 
above method to detect a cyclic attractor whose period is 
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at most p max , but for that purpose, we will define many 
more variables Xi,t, Xi,t,bu ,-,h f° r t € {0,1, ■ ■ ■ , p max — 1}. 

ILP for SIMULTANEOUS ATTRACTOR CONTROL 

Consider that each variable vi has two possibilities, i.e, 

(i) v t is not chosen for a control node (i.e., v, corre- 
sponds to an internal node), 

(ii) Vi is chosen for a control node (i.e., v t becomes 
an external node). 

In order to specify the state of a variable, we introduce 
the following variables and constraints for N\. 
Let Xi be 

x .\Pi if m = 0; 
1 [Zi if Wi = 1 . 

This function can be represented by 

Pi — Wi < Xi <pi + Wi, 
Zi — (1 — Wi) < Xi < Zi + (1 — Wi). 

The above representation means that w, = 1 corre- 
sponds to the case that x t is selected as a control node, 
to which Z[ gives a 0-1 assignment. Similarly, for N 2 of 
the normal cell, another variable s t is defined as 




cfr if Wi = 0; 
Zi if Wi = 1. 



This function can also be represented by 

Cji — Wi < Si < Cji + Wi 
Zi — (1 — Wi) < Si < Zi + (1 — Wi). 

In this representation, we can see that w t = 1 also cor- 
responds to the case that si is selected as a control 
node, and has the same 0-1 assignment. This guarantees 
that these two different BNs are under the same control. 
For this attractor control problem, we assume that the 
score functions for N x and N 2 take the following forms, 
respectively: 

h(x) = • (1 — u>i) ■ Xi 

i 

and 

i 

We can assume without loss of generality that a, > 0 
and ^ < 0. Then, for the first BN, N lt we try to maximize 
the score of the internal nodes, so we need to maximize 
h(x). The score function can be reduced to E, at ■ Ui if 
we add the constraints u t < x t and m, < 1 - w it where u t is 



a binary variable. In terms of maximizing the score func- 
tion g(s) for N 2 , we introduce the additional binary vari- 
able ji such that y t = -/},-. Thus, the score function for g(s) 
becomes • (1 — if,) • (1 — s,). Furthermore, the 
score function can be converted into £, Yi ' r i, if we add 
the constraints r, < (1 - w,) and r, < (1 - s,). 

The original aim is that the minimum score of single- 
ton attractors is maximized for these two networks. 
However, it is difficult to give a direct ILP formalization 
for this problem, because of the Ej-hardness °f tf* e P ro " 
blem, so as in [15], we firstly formalize an ILP to find a 
singleton attractor with the maximum score by choosing 
and controlling m nodes. Incorporating the inequalities 
mentioned above with the ILP formalization for 
ATTRACTOR DETECTION, the following ILP formali- 
zation for SAC is obtained. 

maximize^ oii ■ Ui + y, ■ n 

subject to 

x 'M l -\>(j2 j , 1 ,...j, T s^))-(. k -i). 

for all i € [1 ■ • • n] and b h ■ ■ ■ !>,, s (0, lj'such that f<(b h , ■■■ , b h ) = 1 
-\ - 0, 

for aU i e [1 • • • n] and b t , ■ ■ ■ b it e (0, lj* such that /(ft,, ■ ■ ■ , b it ) = 0 
^E,,,,,,*/".-'.' foraUi€[l...n] 

£E i ,, i ,. V{0 , 1 i>wv 

pi — Wj < Xi <pi + Wi, 

Zj + Wj— 1 < Xi < Zj — Wi + 1, 

Ui <Xi, Ui < 1 - Wi, for each i € [1 ■ ■ ■ n] 

s *, ~\- (E j , 1 ,., t i s( 5 v)) -(*-!)■ 

s '.K\^l(H j , li ...j,^i l ))' 

for all i e [1 • • • n] and b h ■ ■ ■ b h e [0, lj'such that fi(b h , ■■■ , b,J = 1 
Hh,~\ =0. for all i s [1 ■■ n] and b h ■ -b, t s |0, 1}* such that ./}(&*,, ••■ , b,J = 0 
^En^M' foraUi€[l...n] 

f/j — < 5j < C\i + Ifj, 

Zf + If; — 1 < Sj < Zi — Wi + 1, 

Ti < 1 — s„ n < 1 — ifj, for each i € [1 ■ ■ ■ n] 
*i/ Si- Pi- <7i- Zi, w it Ui, ri€{0,l}- for alii € [l--n], 
Xi,\~\ E {0, 1), for all i e [1 ■ ■ ■ n] and b h ■ ■ ■ b k € (0, 1]' 
Hk, ~\ s {0, 1), for all i € [1 ■ ■ n] and !),,-■■ bi, € (0, 1)* 

We denote this ILP formulation as ILP-A. Since sin- 
gleton attractors are not always uniquely determined, 
considering the worst case, it is necessary to maximize 
the minimum score of the singleton attractors. Suppose 
that V = (v^ , ■ ■ ■ , Vi m ) has been selected as the control 
nodes with 0-1 value B' = (b^, ■ ■ ■ , bi m ). These notions 
are also used in the following two problems, so detec- 
tion of the attractor with the minimum score with these 
control nodes can be formulated as follows. Define / = 
{i v i m }. The objective function can be replaced by 
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"Minimize ^a,x; + — 5,)" 

and the constraints replaced m, and r, by 

)Xj = Zi, Si = Zi, u>i = 1, z; = frj for all i e / 
u>i = 0 for all i ^ J. 

The resulting ILP is denoted by ILP-A'. 

From the perspective of avoiding examination of the 
examined (V, B'), we will modify ILP-A. In other words, 
we try to find node-value pairs that are not the same as 
the previously obtained solutions. This can be tackled by 
some additional linear inequalities which ensures that the 
following node-value pairs differ from the previous ones. 
Note that these two networks have the same control, so 
if we can avoid obtaining the previously examined (V, B') 
for N u then we can also get different node-value pairs 
for N 2 . Thus, we shall consider one of the networks, i.e., 
N v Assume that x k = (xf\ xf\ ■■■ , x ( „ ] ) is the kth con- 
trol previously found, where we let x- fe ' = z-^ if w\ = 1, 
and otherwise xf 1 = — 1- Then for each k, we add the fol- 
lowing linear inequality: 

where 5{x, y) is the delta function, 



S{x, y) 



( 1 if x 
jo if i 



= Y 



This inequality guarantees that the following must 
hold for at least one of the control nodes: 

. if x \ k) = 1, either z f +1) = 0 or w f +1] = 0 holds, 
♦ otherwise, either z | fe+1) = l or w| fe+1) = 0 holds. 

We consider the case that one of the control nodes xj^ 
is (i.e., zf ] = 1) for the kth control. If w f +1) = 0 holds for 
the {k + l)th control, then in the next iteration, it will not 
be selected as a control node. Otherwise, w\ k+1 ^ = L so it 
is still a control node. However, the inequality guarantees 
that zf +1 ^ must be 0, which is different from the previous 
value (z| fe ' = 1)- Thus, the above inequality ensures that 
the following obtained set of node-value pairs is different 
from the previous ones. Define ILP-B as this modified 
version in the following context. We can solve the simul- 
taneous attractor control problem for multiple networks 
through the following algorithm. 

1. Repeat steps 2-3. 

2. Find (V, B') yielding the maximum score of a sin- 
gleton attractor using ILP-B where (V, B) is not the 



same as any of the already examined nodes/values 
pairs. "NULL" should be output and halt, if the max- 
imum score is less than 6. 

3. For (V, B'), calculate the singleton attractors with 
the minimum score by ILP-A'. Output (V, B') and 
halt, if the minimum score is no less than 8. 

Note that in the worst case, it may repeat this proce- 
dure exponentially many times, but we expect that the 
procedure will not be repeated so many times, since the 
expected number of singleton attractors (i.e, per (V, B)) 
is small, regardless of the total number of nodes (k). 
How to select the thresholds 8 is an important issue in 
this program. If we know an appropriate threshold in 
advance, we can simply use such a 8. Here, we let 8 be 
1.2m, because the desired attractors almost always exist 
if 8 « 1.2n, and it often occurs in our preliminary com- 
putational experiments that the algorithm cannot find 
the desired attractors if 8 » 1.2k. 

ILP for attractor control under damaging cancer cells 
substantially 

We try to investigate if there exists a control that 
ensures damaging the cancer cell substantially and, 
under this condition, identifying singleton attractors for 
the normal cell. For the ILP-A of ACDC (Attractor 
Control under Damaging Cancer cells substantially), we 
have to make sure that the control damages the cancer 
cell significantly by introducing the following inequality 



J2 YtUt - 



(1) 



A larger ^ signifies that the cancer cell is damaged 
substantially. Here we set the fi as 0.6k. Then, for the 
normal cell, we define the objective function as 

Yi ■ ?i. It should be noted that for this problem it is 

possible that there does not exist any singleton attractor 
for A/2 since the control set must satisfy the inequality 
(1) ensuring significant damage to the cancer cells. In 
terms of ILP-A', we aim at finding the maximum of the 
minimum score of N 2 for the case of having multiple 
singleton attractors. Since we consider the multiple sin- 
gleton attractors for N 2 , we do not need to include con- 
straints about N lt so we can replace r, with 



Si = Zi, u>i = 1, Zi = hi 
u>i = 0 



for all i € I 
for all i g I 



Consequently, the objective function becomes "Mini- 
mize yi^j K J (1 — s i)"- As for ILP-B, to avoid examining 
the previously obtained (V, B') for N 2 , we assume that 
s k = (sf',s 2 k) , . . . ,s { n ] ) is the ^ th control previously 
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found, where we let sf' = s) if w t = 1, otherwise 
s[ fe ' = — 1. Then for each k, we add the 
following linear inequality: 

We can solve this problem using the following algo- 
rithm. 

1. Repeat steps 2-3. 

2. Find (V, B) yielding the maximum score of a sin- 
gleton attractor using ILP-B where (V, B") is not the 
same as any of the already examined nodes/values 
pairs and that the control damages the cancer cell 
significantly. "NULL" should be output and halt, if 
the singleton attractor does not exist. 

3. For (V, B'), calculate the singleton attractors with 
the minimum score by ILP-A'. Output (V, B') and 
halt, if the minimum score is greater than 0 V 

Here, we set 0\ to be 0.6k, because if 0\ « 0.6k, then 
the desired attractor almost always exists, whereas if 
6i » 0.6k, then the desired attractor seldom exists. 

ILP for attractor control under keeping normal cells 
undamaged 

We also try to investigate whether there exists a control 
that ensures not damaging the normal cell substantially 
and under this condition, we find a singleton attractor for 
the cancer cell. In terms of ILP-A for ACKN (Attractor 
Control under Keeping Normal cells undamaged), consid- 
ering this control should not damage the normal cell sub- 
stantially, so we introduce an additional inequality 

I>;n>&. (2 ) 

i 

We set the £2 to be 0.6m, and we define the objective 
function as «i • M,. The singleton attractor for ATi is 

difficult to find for this problem, because the control set 
must satisfy the inequality (2) (not damaging the normal 
cell substantially). Considering that there may exist mul- 
tiple singleton attractors for the cancer cell, and the 
worst case, it is necessary to maximize the minimum 
score of singleton attractors. 

Thus for ILP-A', we do not need to add any con- 
straints about N 2 for ILP-A' so we replace u t by 

Xi = Zi, Wi = 1, z; = hi for all i € I 
u>i = 0 for all i g I 

The objective function is "Minimize ^^ J Q, ! X >- Except 
these modifications, ILP-B is exactly the same as in the 



case of SAC. We can solve this problem through the fol- 
lowing algorithm. 

1. Repeat steps 2-3. 

2. Find (V, B) yielding the maximum score of a sin- 
gleton attractor using ILP-B where (V, B) is not the 
same as any of the already examined nodes/values 
pairs and that the control causes only limited 
damage to the normal cell. "NULL" should be out- 
put and halt, if no singleton attractor for ATi exist. 

3. For (V, BX calculate the singleton attractors with 
the minimum score by ILP-A'. Output ((V, B) and 
halt, if the minimum score is not less than 0 2 . 

For this problem, we set 0 2 as 0.6k because it often 
occurs that the algorithm can find the desired attractor 
if 0 2 « 0.6k, and that it hardly finds the desired attrac- 
tor if 0 2 » 0.6k. 

Results 

Results on simultaneous attractor control 

The proposed method for SIMULTANEOUS ATTRAC- 
TOR CONTROL was applied to randomly generated 
BNs. The following data is according to 10 randomly 
generated pairs of BNs with K = 2, where for each pair 
of BNs, the BNs have the same nodes but different Boo- 
lean functions. 

1. time: average CPU time (seconds) per pair of BNs, 

2. #pos#rep: the number of pairs of BNs where the 
desired attractors can be found, and the average 
number of repeats for each such pair of BNs (recall 
that it is necessary for SAC (Simultaneous Attractor 
Control) to solve ILP instances repeatedly), 

3. #neg#rep: the number of pairs of BNs for which 
the desired attractor does not exist, and the average 
number of repeats for each such pair of BNs. 

We set a, = 1 and J, = 1 for all i G {1, 2, k}, and we 
set the maximum number of repeats to be 20. It is pos- 
sible that this procedure may not decide whether or not 
the desired attractors exist within 20 repeats in some 
cases. The table shows that the number of repeats is 
small even though the number of nodes is large. The 
reason the CPU time for (140, 14) is greater than that 
for (160, 16) is that the latter required fewer repeats. 

Results on attractor control for normal cells under 
damaging cancer cell substantially 

For this problem, we first guarantee that this control 
damages the cancer cell substantially and add the addi- 
tional constraint (1). It is possible that no singleton 
attractor exists for the normal cell under the condition 
that the control set guarantees significant damage to the 
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Table 3 Results on Attractor Control for the Normal Cells 
Given that the Control Damage the Cancer Cells 
Significantly 



n/m 


100/10 


120/12 


140/14 


160/16 


time(sec) 


3.24 


4.47 


6.84 


8.94 


#pos/#rep 


7/5 


6/2.33 


4/3.75 


5/5.8 


#neg/#rep 


1/2 


1/2 


3/2 


3/3 


Table 4 Results on Attractor Control for the Cancer Cells 


Given that the control has limited Damage to the Normal 


Cells 










n/m 


100/10 


120/12 


140/14 


160/16 


time(sec) 


3.12 


4.05 


5.89 


8.21 


#pos/#rep 


8/2.88 


5/2.6 


4/1.5 


4/3.5 


#neg/#rep 


2/2 


4/1 


4/1 


4/1.5 



cancer cells. As shown in Table 3, our proposed method 
is useful also for this problem. 

Results on attractor control for cancer cells under 
keeping normal cells undamaged 

The results also suggest that our proposed approach is 
efficient and effective for solving the problem of attrac- 
tor control of the cancer cell with no damage to the 
normal cell guaranteed. 



In order to verify our proposed method further, we 
applied ILP to some real networks (see Figure 5B) in [16]. 
There are two types of mice in that figure. The left figure is 
for BALB/c, which is more sensitive to radiation and tends 
to become cancer, while the right one is for C57BL/6, 
which is more resistant to radiation and tends not to 
become cancer. However, the nodes for these two networks 
are different and the Boolean functions are not given. Thus, 
we consider utilizing the right figure, which corresponds to 
the normal cell, and further simulating a cancer cell based 
on the normal cell (right figure) to guarantee that both net- 
works have the same nodes (see Figure 2(b)). Furthermore, 
the red and green nodes for the normal cell represent 1 and 
0, respectively. We assign most of the OR nodes to the can- 
cer cell, and we see that the majority of nodes in the single- 
ton attractor of the cancer cell are 1. Similarly, we assign 
most of the AND nodes to the normal cell, and we see that 
the majority of nodes in the singleton attractor of the nor- 
mal cell are 0. This may mean that this is a subnetwork 
activated by cancer. Thus we have obtained the singleton 
attractor for the cancer cell and normal cell. Since our 
objective is to transform the state of the cancer cell into 
that of the normal cell, we try to find a control such that 
most of the 1 nodes in the cancer cell are converted into 0 
nodes. Then we define a score function for both BNs. 
Specifically, we assume % = 1 and a, = 1 for all i, and we set 
n = 36, m = 4 and 0 = 1.2-n. We can see that most of the 
states of the singleton attractor for cancer cells are changed 




(a) normal cell 

Figure 2 Normal cell vs cancer cell 




(b) cancer cell 
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into 0, indicating they are the desired attractors for the can- 
cer cell under this control. Moreover, the CPU time is 0.03 
(sec), which suggests that the proposed method might work 
efficiently for real medium-size networks. 

Discussions 

In this paper, we formulated three novel problems, 
simultaneous attractor control for multiple networks, 
attractor control for normal cells with guaranteed 
damage to cancer cells, and attractor control for cancer 
cells guaranteed not to damage normal cells, and we 
applied an ILP-based approach for solving these pro- 
blems in a unified manner. We further investigated the 
attractor control problem for multiple BNs and vali- 
dated our proposed method for realistic networks. 
Though attractor control problems are Z/>-hard, the 
experimental results have shown the efficiency of our 
proposed method. Furthermore, this method was seen 
to be useful for solving medium-size instances of these 
problems. The method we proposed might not be the 
fastest, but it is easy and simple to implement and, 
furthermore, it has rooms for modifications and exten- 
sions. In particular, the use of non-linear costs for the 
scoring functions is of interest for future work. 
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